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PREFACE 


As  part  of  The  RAND  Corporation's  continuing  study  of  system 
identification,  this  Memorandum  presents  a  method  of  applying  quasi¬ 
linearization  and  digital  computers  to  the  problem  of  simultaneously 
estimating  the  density  of  the  earth's  atmosphere  at  high  altitudes 
and  the  geopotential,  making  use  of  measurements  of  a  satellite's 
motion  through  the  atmosphere.  This  study  will  be  of  interest  to 
those  working  in  the  fields  of  space  science  and  meteorology. 


SUMMARY 


A  new  approach  to  the  problem  of  estimating  the  density  of  the 
earth’s  atmosphere  at  high  altitudes  is  presented.  Using  measurements 
of  a  satellite's  motion  through  the  atmosphere,  we  formulate  the 
problem  as  a  nonlinear  multi-point  boundary  value  problem  which  can 
be  solved  using  quasi linearization  and  high-speed  computers.  A 
general  machine  program  has  been  written,  and  the  results  of  some 
illustrative  preliminary  tests  are  presented.  This  Memorandum  is 
an  introduction  to  a  more  detailed  study. 
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I.  INTRODUCTION 


In  a  recent  paper  King-Hele'*  *  has  described  current  views  of 
the  problem  of  estimating  the  density  of  the  earth's  atmosphere  at 
altitudes  between  200  and  1000  km  employing  measurements  of  earth 
satellites.  In  the  section  on  future  developments,  he  mentions  the 
use  of  measurements  made  during  a  single  revolution  to  gain  improved 
resolution.  The  aim  of  this  Memorandum  is  to  show  how  measurements 
of  displacement  and  velocity  can  be  converted  into  estimates  of  or¬ 
bital  elements,  drag  coefficients  and  potentials.  We  shall  employ  a 

(2) 

quasilinearization  technique.  Earlier  applications  to  orbit  deter¬ 
mination  are  given  in  Refs.  3  and  4.  Our  procedure  takes  full  advan¬ 
tage  of  the  modern  computer's  ability  to  integrate  large  systems  of 
nonlinear  differential  equations  subject  to  multi-point  boundary 
conditions. 


II.  METHOD 


We  shall  present  our  method  through  the  discussion  of  a  simpli¬ 
fied  one-dimensional  example  so  as  not  to  obscure  the  trend  of 
thought  in  a  jungle  of  equations.  Consider  the  motion  of  a  mass  des¬ 
cribed  by  the  differential  equations 

x  =  v  (1) 

v  =  g  -  kv2  .  (2) 

The  parameter  g  is  associated  with  the  gravitational  field,  and  the 
parameter  k  is  a  drag  coefficient.  We  suppose  that  at  various  times 
t^,  i  =  1,  2,  ...,  N,  noisy  observations  are  made  of  the  displacement 
and  velocities,  the  measured  values  at  time  t.  being  b.  and  c.} 

ill 

i  =  1,  2,  ...,  N.  We  wish  to  estimate  the  parameters  g  and  k  and  the 
initial  conditions  (orbital  elements),  x(0)  and  v(0) ,  so  that  we  shall 
minimize  S,  the  sum  of  the  squares  of  the  deviations  between  the  ob¬ 
served  values  and  the  theoretical  values, 


S  = 


(3) 


where  and  \ ^  are  certain  weights.  The  theoretical  values  x(t^) 
and  v(t^)  are  obtained  by  numerically  solving  the  differential  equa¬ 
tions  (?)  and  (2)  with  the  considered  values  of  g  and  k  and  the  con¬ 
sidered  initial  conditions.  The  values  of  these  four  quantities  that 


minimize  S  are  our  estimates. 


We  solve  this  minimization  problem  by  a  numerical  technique  which 

(2) 

has  worked  well  in  many  instances  in  the  past.  '  The  scheme  is  a 

successive  approximation  technique.  We  denote  the  previous  approxi- 

0  0 

mation  of  the  functions  x  and  v  by  x  and  v  ,  for  t^  £  t  ^  t^,  and  the 
new  ones  by  x^  and  v^.  The  previous  estimates  of  g  and  k  are  denoted 
by  and  k^,  and  the  new  ones  by  and  k^.  The  linearized  equations 
for  x^  and  v^  are 


.1  0  ,  ,  1  0.  1 

x  =  v  +(v  -  v  )  =  v 


(4) 


.1  0  ,0/  0.2  ,  ,  „0  0W  1  (h  ,  .  1  (L 

v  =  g  -  k  (v  )  +  (-2k.  v)(v  -  v  )  +  (g  -  g) 


(5) 


.  ,  0.2  n  1  . 0. 

+  (-v  )  (k  -  k  )  . 


These  equations  are  obtained  by  expanding  the  right-hand  sides  of 
Eqs.  (1)  and  (2)  about  x  =  x^,  v  =  v^,  k  =  k^  and  g  =  g^  and  retaining 
only  the  terms  that  are  linear  in  x  ,  v  ,  k^  and  g^.  We  wish  to  deter¬ 
mine  the  functions  x^  and  v^  and  the  constants  g^  and  k^  so  as  to  mini¬ 
mize  the  sum 


This  problem  may  be  solved  computationally  at  the  expense  of  solving 
a  few  initial  value  problems  by  the  method  of  complementary  functions, 
a  task  for  which  modern  computers  are  admirably  suited. 

First  we  produce  numerically  the  solutions  of  the  homogeneous 
equat ions 
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hll  h21’ 


K,  =  (-2k  vu)  h01,  t.  <;  t  £  tM, 


hll(V  1 5 


h21(tx)  =  0, 


^12  h22 


b22  =  ^"2k  V  ^  h22s  *  fc  *  CN> 


h12^1^  =  0> 


h22Ctl)  =  U 


Then  we  numerically  produce  the  particular  solutions  p^  and  p^ 


Pi  =  P2> 

.  0  .0,  0  2  ...  0  0,  ,  0.  0  ,  .  0N  2  ,  ,0.  / t  c\ 

P0  =  8  -  k  v.v  )  -  (2k  v  )  (p,  -  v  )  -  g  +  (-v  )  (-k  )  ,  (16) 


on  the  interval  t.  ^  t  s  t,T  using  the  initial  conditions 

1  N 


p^tp  =  o, 


P2 j_)  —  0  . 
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And  finally  we  produce  numerically  the  particular  solutions 


qx  -  q2> 

q1(t1)  “  o, 

(19) 

32  =  (-2k°v°)  q2  +  1, 

q2(t1)  =  °> 

(20; 

and 

rl  =  r2 5 

H-* 

/■N 

r r 

h-4 

v-/ 

II 

O 

(21) 

/  „0  0,  .  ,  0n2 
r2  =  (-2k  v  )  r2  +  (-v  )  , 

r2(tx)  =  0, 

all 

for  t  £  t  ^  t2« 

The  general  solution  to  Eqs.  (4)  and 

(5)  may  now  be 

written  in 

the 

form 

x1  =  ck1  hn  -h  3  h12  +  p^  + 

g  qx  +  k  rx 

(22) 

V1  -  a  h2l  +  f>  h22  +  p2  + 

g  q2  +  k  r2  . 

(23) 

The  constants  O'  and  {3  are  the  new  estimates  of  the  orbital  elements: 


x1(t1)  =  a 

(24) 

v1(t1)  =  3  . 

(25) 

To  find  the  values  of  the  parameters  o',  {3,  g  and  k1  which  mini¬ 
mize  the  sum  of  the  squares  of  the  deviations  we  return  to  Eq.  (6) 
and  replace  x^(t^)  and  v^(t^)  by  their  equivalents  from  Eqs.  (22)  and 
(23).  In  this  way  we  see  that  S  is  quadratic  in  the  parameters  a,  {3, 
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g1  and  k1.  We  find  the  values  of  these  parameters  which  minimize 
by  solving  the  linear  system  of  algebraic  equations 


(26) 


(27) 

(28) 

(29) 


Once  the  minimizing  values  of  a,  3,  g  and  k  have  been  determined, 
we  return  to  Eqs.  (22)  and  (23)  and  produce  the  numerical  values  of 
the  functions  x^  and  v^  on  the  interval  t^  £  t  £  t^.  This  completes 
a  basic  cycle  of  the  iterative  calculation.  The  rapid  quadratic  con¬ 
vergence  which  can  be  expected  is  discussed  in  Ref.  2. 

The  initial  approximations  for  the  constants  g  and  k  and  the 
functions  x  and  v  (for  £  t  £  t^)  may  be  obtained,  e.g. ,  by  using 
the  best  estimates  of  g  and  k  available  and  then  integrating  Eqs.  (1) 
and  (2)  from  t^  to  t^,  using  the  best  available  estimates  of  x(0) 
and  v(0). 
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III.  SOME  NUMERICAL  EXPERIMENTS 

First  we  present  the  results  of  an  experiment  designed  to  show 
the  rapidity  of  convergence  of  the  technique.  We  integrated  Eqs.  (1) 
and  (2)  using 

g  -  32.2  (30) 

k  =  32.2  X  10"3  =  0.0322  (31) 

and 


x(0)  =0.0 

(32) 

«•> 

o 

• 

o 

II 

o 

> 

(33) 

the  integration  extending  from  t  =  0  to  t  =  1.0,  The  results  are 
given  in  columns  3  and  5  of  Table  1,  and  provide  us  with  a  set  of  dis¬ 
placement  and  a  set  of  velocity  measurements  for  t  =  .1,  .2,  ...,  1.0. 
The  estimation  problem  consists  of  assuming  these  twenty  observations 
as  given  and  seeking  to  uncover  the  values  of  the  parameters  g  and  k 
as  well  as  the  initial  conditions  x(0)  and  v(0). 

We  employed  the  method  sketched  in  Section  II,  assuming  our 
initial  estimates  of  the  unknowns  to  be 


x°(0)  =  0.1 

(34) 

v°(0)  =  -0.1 

(35) 

g°  =  30.0 

(36) 

k°  =  0.003  . 


(37) 


8 


Table  1 

OBSERVED  AND  THEORETICAL  VALUES  OF 


x  Observations 


Time 

X 

b. 

l 

0.1 

0. 16072253E 

00 

0. 16072254E 

00 

0.2 

0. 63959687E 

00 

0. 63959709E  00 

0.3 

0. 14270092E 

01 

0. 14270096E 

01 

0.4 

0. 25077758E 

01 

0. 25077762E 

01 

0.5 

0. 38622596E 

01 

0. 38622599E 

01 

0.6 

0. 54676883E 

01 

0. 54676887E 

01 

0.7 

0. 72995026E 

01 

0.  72995033E 

01 

0.8 

0,  93325923E 

01 

0. 93325937E  01 

0.9 

0.11542327E 

02 

0.  11542329E 

02 

1.0 

0. 13905332E 

02 

0. 13905336E 

02 

The  E  notation  refers  to  exponents 
=  3.2089171. 


-V 


DISPLACEMENTS  AND  VELOCITIES 


v 

0. 32089146E  01 
0. 63524209E  01 
0. 93703303E  01 
0. 12212048E  02 
0. 14839446E  02 
0. 17227990E  02 
0. 19366275E  02 
0.2 125434 IE  02 
0.22901304E  02 
0.24322771E  02 


v  Observations 
c . 

_ x _ 

0. 32089171E  01 
0. 63524224E  01 
0. 93703309E  01 
0.  122 12 049E  02 
0. 14839446E  02 
0. 17227991E  02 
0. 19366278E  02 
0. 21254349E  02 
0.22901317E  02 
0. 24322791E  02 


of  10,  i.e.  ,  0.32089171E  01 
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Notice  that  the  initial  approximation  of  k  is  incorrect  by  a  factor  greater 
than  ten.  Our  initial  estimates  of  the  functions  x(t)  and  v(t)  were 
obtained  by  integrating  Eqs.  (1)  end  (2)  subject  to  the  conditions  in 
Eqs.  (34)  through  (37).  results  are  displayed  in  Table  2.  They 

show  that  excellent  estimates  are  obtained  by  the  '’^nrnximation. 

The  second  and  fourth  columns  of  Table  1  contain  the  values  of  the  dis¬ 
placements  and  velocities  obtained  using  those  estimates  in  Eqs.  (1) 
and  (2 ) . 

Now  we  turn  to  the  results  of  an  experiment  involving  noisy  ob¬ 
servations.  The  values  of  the  velocities  used  in  the  earlier  experi¬ 
ment  were  corrupted  by  adding  about  two  percent  gaussian  noise  to  each 
observation.  These  new  noisy  observations  are  listed  in  column  5  of 
Table  3.  The  displacement  observations  are  the  same  as  before,  and  the 
value  of  g  is  kept  at  32. 2.  We  found  the  minimizing  values  of  the  un¬ 
knowns  to  be 

x(0)  =  0.928198  X  10_1  (38) 

v(0)  =  -0.270570  (39) 

k  =  0.298726  X  10"1  .  (40) 

This  serves  as  a  warning  by  showing  that  quite  small  errors  in  the  ob¬ 
servations  can  lead  to  significantly  larger  errors  in  the  estimates. 

The  value  of  using  numerical  experiments  such  as  these  to  aid  in 
estimating  the  accuracy  and  number  of  the  measurements  required  to 
obtain  the  desired  accuracy  in  the  estimates  is  evident. 

Results  such  as  these  are  obtained  in  an  execution  time  of  a  few 
seconds  on  an  IBM  7044. 


Table  3 


THE  FIT  TO  THE  NOISY  OBSERVATIONS 


Time 

X 

x  Observations 
b. 

i 

V 

v  Observations 

c . 

l 

0.1 

0.226581E 

00 

0. 160722E 

00 

0.294152E 

01 

0.320540E 

01 

0.2 

0. 679259E 

00 

0. 639597E 

00 

0. 609782E 

01 

0. 627940E 

01 

0.3 

0. 144238E 

01 

0. 142701E 

01 

0. 914155E 

01 

0. 882816E 

01 

0.4 

0. 250218E 

01 

0.250778E 

01 

0. 120240E 

02 

0.  121049E 

02 

0.5 

0.3 84052 E 

01 

0.  38  622  6E 

01 

0. 147071E 

02 

0. 144243E 

02 

0.6 

0. 543608E 

01 

0. 546769E 

01 

0. 171651E 

02 

0. 170832E 

02 

0.7 

0. 726556E 

01 

0. 729950E 

01 

0.  193840E 

02 

0. 192415E 

02 

0.8 

0.  930481E 

01 

0. 933259E 

01 

0.213608E 

02 

0. 212337E 

02 

0.9 

0.  115298E 

02 

0. 115423E 

02 

0.  231011E 

02 

0. 236843E 

02 

1.0 

0.  139176E 

02 

0. 139053E 

02 

0. 246175E 

02 

0. 243658E 

02 
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IV.  DISCUSSION 

A  general  program  in  FORTRAN  IV  has  been  written  which  solves 
nonlinear  multi-point  boundary  value  problems  using  the  method  sketched 
in  Section  II.  It  is  available  from  the  authors.  We  plan  to  investi¬ 
gate  its  utility  in  estimating  parameters  of  the  geopotential  and  drag 
coefficient  from  satellite  measurements. 

Previous  experience  plus  experiments  similar  to  those  reported 
here  show  that  small  parameters  in  differential  equations  can  be  esti¬ 
mated  successfully  using  experimental  observations  on  the  solutions. 

It  must  not  be  thought,  however,  that  matters  are  completely 
routine.  Various  problems  such  as  instability  of  the  linearized 
equations,  limited  memory,  and  ill-conditioning  of  the  linear  algebraic 
equations  can  and  do  arise.  Methods  for  overcoming  these  difficulties 
are  given  in  Refs.  2  and  5. 
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